Fit condition model with environmental and biological predictors

Fit main model (see exploratory scripts and model comparison), visualize results.

Load packages

library(tidyverse)
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(RColorBrewer)
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2) 
library(ggcorrplot)
library(ggridges)
library(brms)
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) 
library(sdmTMBextra) 

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")

Code for map plots

# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
sf::sf_use_s2(FALSE)
# https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data

swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

# Define plotting theme for main plot
theme_plot <- function(base_size = 11, base_family = "") {
  theme_light(base_size = base_size, base_family = "") +
    theme(
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 9, colour = 'gray10', margin = margin(b = 1, t = 1)),
      strip.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 11, base_family = "") {
  theme_light(base_size = base_size, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 7),
        strip.text = element_text(size = 8, colour = 'gray10', margin = margin(b = 1, t = 1)),
        strip.background = element_rect(fill = "gray95"),
        legend.direction = "horizontal",
        legend.margin = margin(1, 1, 1, 1),
        legend.box.margin = margin(0, 0, 0, 0),
        legend.key.height = unit(0.4, "line"),
        legend.key.width = unit(2, "line"),
        legend.spacing.x = unit(0.1, 'cm'),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())
}

# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2

plot_map_raster <- 
ggplot(swe_coast_proj) + 
  xlim(xmin2, xmax2) +
  ylim(ymin2, ymax2) +
  labs(x = "Longitude", y = "Latitude") +
  geom_sf(size = 0.3) + 
  theme_plot() +
  guides(colour = guide_colorbar(title.position = "top", title.hjust = 0.5),
         fill = guide_colorbar(title.position = "top", title.hjust = 0.5))

Read data

# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(2_2).csv")

d <- bind_rows(d1, d2)

unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE

# Calculate standardized variables
d <- d %>% 
  mutate(ln_length_cm = log(length_cm),
         ln_weight_g = log(weight_g),
         year_ct = year - mean(year),
         biomass_her_sc = biomass_her,
         biomass_her_sd_sc = biomass_her_sd,
         biomass_spr_sc = biomass_spr,
         biomass_spr_sd_sc = biomass_spr_sd,
         density_cod_sc = density_cod,
         density_cod_rec_sc = density_cod_rec,
         density_fle_sc = density_fle,
         density_fle_rec_sc = density_fle_rec,
         density_saduria_sc = density_saduria,
         density_saduria_rec_sc = density_saduria_rec,
         depth_sc = depth,
         depth_rec_sc = depth_rec,
         oxy_sc = oxy,
         oxy_rec_sc = oxy_rec,
         temp_sc = temp,
         temp_rec_sc = temp_rec,
         year_f = as.factor(year))  %>%
  mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
              "density_cod_sc", "density_cod_rec_sc", 
              "density_fle_sc", "density_fle_rec_sc", 
              "density_saduria_sc", "density_saduria_rec_sc", 
              "depth_sc", "depth_rec_sc", 
              "oxy_sc", "oxy_rec_sc", 
              "temp_sc", "temp_rec_sc"
              ),
            ~(scale(.) %>% as.vector)) %>% 
  mutate(year = as.integer(year))

unique(is.na(d))
#>      weight_g length_cm  year   lat   lon ices_rect sub_div depth   oxy  temp
#> [1,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [2,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [3,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [4,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [5,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [6,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [7,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#>          X     Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [2,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [3,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [4,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [5,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [6,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [7,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#>      density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,]           FALSE               FALSE           FALSE           FALSE
#> [2,]           FALSE               FALSE           FALSE           FALSE
#> [3,]           FALSE               FALSE           FALSE           FALSE
#> [4,]           FALSE               FALSE           FALSE           FALSE
#> [5,]           FALSE                TRUE            TRUE            TRUE
#> [6,]           FALSE                TRUE            TRUE            TRUE
#> [7,]           FALSE               FALSE           FALSE           FALSE
#>      depth_rec oxy_rec temp_rec biomass_spr biomass_her biomass_spr_sd
#> [1,]     FALSE   FALSE    FALSE       FALSE       FALSE          FALSE
#> [2,]     FALSE   FALSE    FALSE        TRUE        TRUE          FALSE
#> [3,]     FALSE   FALSE    FALSE        TRUE        TRUE           TRUE
#> [4,]     FALSE   FALSE    FALSE       FALSE       FALSE          FALSE
#> [5,]      TRUE    TRUE     TRUE        TRUE        TRUE          FALSE
#> [6,]      TRUE    TRUE     TRUE        TRUE        TRUE          FALSE
#> [7,]     FALSE   FALSE    FALSE        TRUE        TRUE          FALSE
#>      biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [2,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [3,]           TRUE        FALSE       FALSE   FALSE           TRUE
#> [4,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [5,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [6,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [7,]          FALSE        FALSE       FALSE   FALSE           TRUE
#>      biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,]             FALSE          FALSE             FALSE          FALSE
#> [2,]             FALSE           TRUE             FALSE          FALSE
#> [3,]              TRUE           TRUE              TRUE          FALSE
#> [4,]             FALSE          FALSE             FALSE          FALSE
#> [5,]             FALSE           TRUE             FALSE          FALSE
#> [6,]             FALSE           TRUE             FALSE          FALSE
#> [7,]             FALSE           TRUE             FALSE          FALSE
#>      density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,]              FALSE          FALSE              FALSE              FALSE
#> [2,]              FALSE          FALSE              FALSE              FALSE
#> [3,]              FALSE          FALSE              FALSE              FALSE
#> [4,]              FALSE          FALSE              FALSE              FALSE
#> [5,]               TRUE          FALSE               TRUE              FALSE
#> [6,]               TRUE          FALSE               TRUE              FALSE
#> [7,]              FALSE          FALSE              FALSE              FALSE
#>      density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc oxy_rec_sc temp_sc
#> [1,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [2,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [3,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [4,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#> [5,]                   TRUE    FALSE         TRUE  FALSE       TRUE   FALSE
#> [6,]                   TRUE    FALSE         TRUE  FALSE       TRUE   FALSE
#> [7,]                  FALSE    FALSE        FALSE  FALSE      FALSE   FALSE
#>      temp_rec_sc year_f
#> [1,]       FALSE  FALSE
#> [2,]       FALSE  FALSE
#> [3,]       FALSE  FALSE
#> [4,]       FALSE  FALSE
#> [5,]        TRUE  FALSE
#> [6,]        TRUE  FALSE
#> [7,]       FALSE  FALSE

unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE
unique(is.na(d$density_cod))
#> [1] FALSE

# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
                   density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
                   density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
                   oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)

# Sample size
nrow(d)
#> [1] 94295

# stdev of oxygen
sd(d$oxy)
#> [1] 1.805214


# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
                             "biomass_spr_sc", "biomass_spr_sd_sc",
                             "density_cod_sc", "density_cod_rec_sc", 
                             "density_fle_sc", "density_fle_rec_sc", 
                             "density_saduria_sc", "density_saduria_rec_sc", 
                             "depth_sc", "depth_rec_sc",
                             "oxy_sc", "oxy_rec_sc",
                             "temp_sc", "temp_rec_sc")

corr <- round(cor(d_cor), 1)

ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.3))

ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)

# Plot data
d_plot <- d %>% 
  group_by(year, X, Y) %>% 
  summarise(n = n())

plot_map_raster + 
  geom_point(data = d_plot, aes(X*1000, Y*1000, size = n, color = n), alpha = 0.5) + 
  scale_color_viridis() + 
  scale_size(range = c(.1, 3), name = "N per haul") + 
  facet_wrap(~year, ncol = 5) + 
  guides(size = "none") + 
  theme_facet_map()

ggsave("figures/supp/condition/spatial_cond_data.png", width = 6.5, height = 8.5, dpi = 600)

Read the prediction grids

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")

pred_grid <- bind_rows(pred_grid1, pred_grid2)

unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE

pred_grid <- pred_grid %>% mutate(year = as.integer(year),
                                  year_f = as.factor(year))

# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec,
                oxy, oxy_rec, 
                temp, temp_rec) %>%
  mutate_all(~mean(.)) %>%
  distinct(.keep_all = TRUE)

data_stdev <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec, 
                oxy, oxy_rec, 
                temp, temp_rec
                ) %>%
  mutate_all(~sd(.)) %>%
  distinct(.keep_all = TRUE)

# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles in the sub division.
# Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so 
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>% 
  group_by(year) %>% 
  mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
         median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
  mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
         biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>% 
  group_by(year, sub_div) %>% 
  mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
         median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>% 
  # Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles 
  # with the sub-division value divided by the number of rectangles in that sub division
  group_by(year, sub_div) %>% 
  mutate(n_rect = length(unique(ices_rect))) %>% 
  ungroup() %>% 
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))

pred_grid <- pred_grid %>%
  mutate(ln_length_cm = log(1),
         year_ct = year - mean(year),
         biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
         biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
         biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
         biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
         density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
         density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
         density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
         density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
         density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
         density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
         depth_sc = (depth - data_means$depth) / data_stdev$depth,
         depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
         oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
         oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
         temp_sc = (temp - data_means$temp) / data_stdev$temp,
         temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
         )

# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na() %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na(biomass_spr_sd) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

unique(is.na(pred_grid))
#>          X     Y depth  year    id oxy_q3   oxy temp_q3  temp   lon   lat
#> [1,] FALSE FALSE FALSE FALSE FALSE  FALSE FALSE   FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE  FALSE FALSE   FALSE FALSE FALSE FALSE
#>      sub_div ices_rect density_saduria biomass_spr biomass_her biomass_spr_sd
#> [1,]   FALSE     FALSE           FALSE       FALSE       FALSE          FALSE
#> [2,]   FALSE     FALSE           FALSE       FALSE       FALSE          FALSE
#>      biomass_her_sd density_cod density_fle depth_rec temp_rec oxy_rec
#> [1,]          FALSE       FALSE       FALSE     FALSE    FALSE   FALSE
#> [2,]          FALSE       FALSE       FALSE     FALSE    FALSE   FALSE
#>      density_cod_rec density_fle_rec density_saduria_rec year_f
#> [1,]           FALSE           FALSE               FALSE  FALSE
#> [2,]           FALSE           FALSE               FALSE  FALSE
#>      median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,]                               FALSE                                 FALSE
#> [2,]                               FALSE                                 FALSE
#>      median_biomass_sprat_across_rect_in_sub_div
#> [1,]                                       FALSE
#> [2,]                                        TRUE
#>      median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,]                                         FALSE  FALSE        FALSE   FALSE
#> [2,]                                          TRUE  FALSE        FALSE   FALSE
#>      biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,]          FALSE             FALSE          FALSE             FALSE
#> [2,]          FALSE             FALSE          FALSE             FALSE
#>      density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,]          FALSE              FALSE          FALSE              FALSE
#> [2,]          FALSE              FALSE          FALSE              FALSE
#>      density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#> [2,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#>      oxy_rec_sc temp_sc temp_rec_sc
#> [1,]      FALSE   FALSE       FALSE
#> [2,]      FALSE   FALSE       FALSE

pred_grid %>% 
  drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

Make spde mesh

spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 100,  
                  type = "kmeans", seed = 42)

# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

ptm <- proc.time()
mnull <- sdmTMB(
  formula = ln_weight_g ~ ln_length_cm,
  data = d,
  time = NULL,
  mesh = spde,
  family = student(link = "identity", df = 5),
  spatiotemporal = "off",
  spatial = "off",
  silent = TRUE,
  reml = FALSE
  )
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.5.4
#> Current Matrix version is 1.5.3
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
#> glmmTMB was built with TMB version 1.9.3
#> Current TMB version is 1.9.4
#> Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
proc.time() - ptm
#>    user  system elapsed 
#>   1.971   0.311   2.350

sanity(mnull)
#> ✖ Non-linear minimizer did not converge: do not trust this model!
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Extract average a and b across the total dataset
a <- as.numeric(tidy(mnull)[1, "estimate"])
b <- as.numeric(tidy(mnull)[2, "estimate"])
a
#> [1] -4.638609
b
#> [1] 2.984044

d$weight_g_avg_pred <- exp(a)*d$length_cm^b

plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")

d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)

Fit main model

ptm <- proc.time()
mfull <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
    biomass_spr_sc + biomass_spr_sd_sc +
    density_cod_sc + density_cod_rec_sc + 
    density_fle_sc + density_fle_rec_sc + 
    density_saduria_sc + density_saduria_rec_sc + 
    depth_sc + depth_rec_sc + 
    oxy_sc + oxy_rec_sc + 
    temp_sc + temp_rec_sc,
    data = d,
  time = "year",
  mesh = spde, 
  family = student(link = "identity", df = 5),
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1)#,
  # priors = sdmTMBpriors(b = normal(c(rep(NA, 27), rep(0, 16)),
  #                                  c(rep(NA, 27), rep(1, 16))))
  )
system("say Just finished!")
proc.time() - ptm
#>    user  system elapsed 
#> 208.426   8.873 221.698

sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))

tidy(mfull, effects = "ran_par", conf.int = TRUE) 
#> Standard errors intentionally omitted because they have been calculated in log
#> space.

# Arrange coefficients
tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(desc(estimate))
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.5.4
#> Current Matrix version is 1.5.3
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> filter: removed 27 rows (63%), 16 rows remaining


tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(estimate)
#> filter: removed 27 rows (63%), 16 rows remaining


# Average magnitude of fixed effects?
tidy(mfull, conf.int = TRUE) %>%
  filter(!grepl("year", term)) %>%
  mutate(abs_coef = abs(estimate)) %>% 
  summarise(mean_abs_coef = mean(abs_coef))
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: new variable 'abs_coef' (double) with 16 unique values and 0% NA
#> summarise: now one row and one column, ungrouped


tidy(mfull, conf.int = TRUE, effects = "ran_par") %>%
  filter(term %in% c("sigma_O", "sigma_E")) %>%
  summarise(mean_estimate = mean(estimate))
#> Standard errors intentionally omitted because they have been calculated in log
#> space.
#> filter: removed 2 rows (50%), 2 rows remaining
#> summarise: now one row and one column, ungrouped

# This model is for the sensitivity analysis, comparing estimates between different model subsets etc. These tend to fit better with only a spatiotemporal term (spatial stdev tends to blow up), so I want them to be as comparable as possible
ptm <- proc.time()
mfull_spatiotemporal <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
    biomass_spr_sc + biomass_spr_sd_sc +
    density_cod_sc + density_cod_rec_sc + 
    density_fle_sc + density_fle_rec_sc + 
    density_saduria_sc + density_saduria_rec_sc + 
    depth_sc + depth_rec_sc + 
    oxy_sc + oxy_rec_sc + 
    temp_sc + temp_rec_sc,
    data = d,
  time = "year",
  mesh = spde, 
  family = student(link = "identity", df = 5),
  spatiotemporal = "iid",
  spatial = "off",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1)
  )
system("say Just finished!")
proc.time() - ptm
#>    user  system elapsed 
#> 148.438   6.709 157.957

sanity(mfull_spatiotemporal)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

tidy(mfull_spatiotemporal, conf.int = TRUE) %>% arrange(desc(abs(estimate)))

tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE) 
#> Standard errors intentionally omitted because they have been calculated in log
#> space.



# Extract random and fixed coefficients from the full model
mfull_est_st <- bind_rows(tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE) %>%
                            filter(term == "sigma_E"),
                          tidy(mfull_spatiotemporal, effects = "fixed", conf.int = TRUE)  %>% 
                            filter(!grepl('year', term))) %>%
  
  mutate(term = factor(term))
#> Standard errors intentionally omitted because they have been calculated in log
#> space.
#> filter: removed 2 rows (67%), one row remaining
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: converted 'term' from character to factor (0 new NA)

# Extract coefficients and save as csv
mfull_est_st <- mfull_est_st %>% 
  mutate(group = "Herring",
         group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
         group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
         group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
         group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
         group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
         group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
         group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
         group = ifelse(term == "sigma_E", "Random", group),
         ) %>% 
  mutate(scale = "Subdivision",
         scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
                                    "density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
         scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
                                    "depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
         scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
         ) %>% 
  mutate(term2 = ifelse(term == "sigma_E", 2, 1))
#> mutate: new variable 'group' (character) with 9 unique values and 0% NA
#> mutate: new variable 'scale' (character) with 4 unique values and 0% NA
#> mutate: new variable 'term2' (double) with 2 unique values and 0% NA

# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est_st, "output/mfull_est_st.csv")

# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mfull, "output/mfull.rds")

samps <- sdmTMBextra::predict_mle_mcmc(mfull, mcmc_iter = 201, mcmc_warmup = 200)

mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_samples = samps)

png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()

d$residuals <- mcmc_res

Calculate approximate \(R^2\) manually using a modifed `r2.sdmTMB`

b_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))

VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
 
# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)
X2 <- mfull$tmb_data$X_ij

first_index <- match('year_f2019', names(as.data.frame(mfull$tmb_data$X_ij[[1]]))) + 1
last_index <- ncol(mfull$tmb_data$X_ij[[1]])

X2 <- matrix(unlist(X2), ncol = length(b_fe))[, first_index:last_index] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects

b_ran <- tidy(mfull, "ran_par")
#> Standard errors intentionally omitted because they have been calculated in log
#> space. Confidence intervals are available with `conf.int = TRUE`.
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
 
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)

VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance

# Marginal R2s
denominator <- VarF + VarO + VarE + VarR
denominator_no_yr <- VarF2 + VarO + VarE + VarR

marg <- VarF/denominator
marg_no_yr <- VarF2/denominator_no_yr

out <- data.frame(
    marginal = marg,
    partial_spatial = VarO/denominator,
    partial_spatiotemporal = VarE/denominator,
    conditional = (VarF + VarO + VarE)/denominator,
    all_random = (VarO + VarE)/denominator,
    marginal_random_ratio = marg / ((VarO + VarE)/denominator),
    random_marginal_ratio = ((VarO + VarE)/denominator) / marg
  )
out


out_no_yr <- data.frame(
    marginal = marg_no_yr,
    partial_spatial = VarO/denominator_no_yr,
    partial_spatiotemporal = VarE/denominator_no_yr,
    conditional = (VarF2 + VarO + VarE)/denominator_no_yr,
    all_random = (VarO + VarE)/denominator_no_yr,
    marginal_random_ratio = marg_no_yr / ((VarO + VarE)/denominator_no_yr),
    random_marginal_ratio = ((VarO + VarE)/denominator_no_yr) / marg_no_yr
  )
out_no_yr

Plot residuals

# Residuals on map
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals), size = 0.5) +
  scale_colour_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  labs(color = "residuals") +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 8.5, dpi = 600)

# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
  geom_point(size = 0.1, alpha = 1) +
  geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
  labs(x = "length [cm]", y = "Residuals", color = "") +
  stat_smooth(color = "black", size = 0.5) +
  facet_wrap(~year, ncol = 5) +
  scale_color_gradient2() +
  theme_plot()
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.


ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 8.5, dpi = 600)

ggplot(d, aes(year, length_cm, z = residuals)) +
  stat_summary_2d(bins = 40) +
  scale_fill_gradient2()


# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)

Annual index using `get_index_sims`

pred_avg_sim <- predict(mfull, newdata = pred_grid %>% filter(depth > 10 & depth < 110), nsim = 500)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining

# Weighted version. We could do it directly in the get_index_sims or manually in the prediction grid (but then we don't get an CIs). This is so that we can divide by the sum of weights (needs to be done by year)

weight_sum <- pred_grid %>% 
  filter(depth > 10 & depth < 110) %>% 
  group_by(year) %>% 
  summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped

index_avg_sim <- get_index_sims(pred_avg_sim,
                                area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>% 
  left_join(weight_sum) %>% 
  mutate(
  est = est / density_cod,
  lwr = lwr / density_cod,
  upr = upr / density_cod
  )
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining with `by = join_by(year)`
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)

index_avg_sims <- get_index_sims(pred_avg_sim,
                                 area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
                                 return_sims = TRUE) %>% 
  left_join(weight_sum) %>% 
  mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining with `by = join_by(year)`
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA

ggplot(index_avg_sim, aes(year, est, ymin = lwr, ymax = upr)) +
  geom_ribbon(alpha = 0.2, color = NA) +
  geom_line(size = 1) + 
  geom_line(data = filter(index_avg_sims, .iteration < 26),
            aes(year, est, group = .iteration), inherit.aes = FALSE, alpha = 0.3)
#> filter: removed 12,825 rows (95%), 675 rows remaining


# Make the un-weighted index
ncells <- filter(pred_grid, year == max(pred_grid$year) & depth > 10 & depth < 110) %>% nrow()
#> filter: removed 250,890 rows (97%), 7,689 rows remaining

index_avg_sim_uw <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

index_avg_sims_uw <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
                                    return_sims = TRUE) # This is for calculation on samples!
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

Get sims from prediction with only fixed effects (-year) and random effects

pred_grid_fe <- pred_grid %>% 
  mutate(year_f = as.factor(1993))
#> mutate: changed 249,002 values (96%) of 'year_f' (0 new NA)
  
pred_avg_sim_fe <- predict(mfull, newdata = pred_grid_fe %>% filter(depth > 10 & depth < 110), nsim = 500, re_form = NA)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining

weight_sum <- pred_grid %>% 
  filter(depth > 10 & depth < 110) %>% 
  group_by(year) %>% 
  summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped

index_avg_sim_fe <- get_index_sims(pred_avg_sim_fe,
                                   area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>% 
  left_join(weight_sum) %>% 
  mutate(
  est = est / density_cod,
  lwr = lwr / density_cod,
  upr = upr / density_cod
  )
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining with `by = join_by(year)`
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)

index_avg_sims_fe <- get_index_sims(pred_avg_sim_fe,
                                    area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
                                    return_sims = TRUE) %>% 
  left_join(weight_sum) %>% 
  mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining with `by = join_by(year)`
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA

# With random effects but no fixed year effect
# pred_grid_re <- pred_grid %>% 
#   mutate(year_f = as.factor(1993))
#   
# pred_avg_sim_re <- predict(mfull, newdata = pred_grid_re %>% filter(depth > 10 & depth < 110), nsim = 500)
# 
# weight_sum <- pred_grid %>% 
#   filter(depth > 10 & depth < 110) %>% 
#   group_by(year) %>% 
#   summarise(density_cod = sum(density_cod))
# 
# index_avg_sim_re <- get_index_sims(pred_avg_sim_re,
#                                    area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>% 
#   left_join(weight_sum) %>% 
#   mutate(
#   est = est / density_cod,
#   lwr = lwr / density_cod,
#   upr = upr / density_cod
#   )
# 
# index_avg_sims_re <- get_index_sims(pred_avg_sim_re,
#                                     area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
#                                     return_sims = TRUE) %>% 
#   left_join(weight_sum) %>% 
#   mutate(est = .value / density_cod)

# Plot all together
index_comp_sim <- bind_rows(#index_avg_sim_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
                            index_avg_sim_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
                            index_avg_sim %>% mutate(Prediction = "Full"))
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA

index_comp_sims <- bind_rows(#index_avg_sims_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
                             index_avg_sims_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
                             index_avg_sims %>% mutate(Prediction = "Full")) %>%
  filter(.iteration < 26)
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> filter: removed 25,650 rows (95%), 1,350 rows remaining

ggplot(index_comp_sim, aes(year, est, ymin = lwr, ymax = upr, color = Prediction, fill = Prediction)) +
  geom_ribbon(alpha = 0.2, color = NA) +
  geom_line(size = 1) + 
  geom_line(data = filter(index_comp_sims, .iteration < 26),
            aes(year, est, group = interaction(.iteration, Prediction), color = Prediction), inherit.aes = FALSE, alpha = 0.3) +
  NULL
#> filter: no rows removed

Difference in Le Cren condition factor between the end and start of time period using sims

# Reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>% 
  dplyr::select(-.value, -density_cod) %>% 
  filter(year %in% c(1993, 2019)) %>% 
  pivot_wider(names_from = year, values_from = est) %>% 
  mutate("2019-1993" = `2019` - `1993`) %>% 
  pivot_longer(2:4) %>% 
  ungroup()
#> filter: removed 12,500 rows (93%), 1,000 rows remaining
#> pivot_wider: reorganized (year, est) into (1993, 2019) [was 1000x3, now 500x3]
#> mutate: new variable '2019-1993' (double) with 500 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 500x4, now 1500x3]
#> ungroup: no grouping variables

# Plot these
p1 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, !name == "2019-1993"),
               aes(value, fill = name), alpha = 0.6) +
  scale_fill_brewer(palette = "Dark2") +
  coord_cartesian(expand = 0)
#> filter: removed 500 rows (33%), 1,000 rows remaining

p2 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, name == "2019-1993"),
               aes(value), alpha = 0.6, fill = "grey") +
  geom_vline(xintercept = 0, linetype = 2) 
#> filter: removed 1,000 rows (67%), 500 rows remaining

p1 / p2


# Calculate quantiles for each level (to go in the results section)
year_diff %>% 
  group_by(name) %>% 
  summarise(median = quantile(value, probs = 0.5),
            upr97.5 = quantile(value, probs = 0.975),
            lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped


# Percent change in condition factor
year_diff %>% 
  group_by(name) %>% 
  pivot_wider(names_from = name, values_from = value) %>%
  mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>% 
  ungroup() %>% 
  summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
            median = quantile(perc_change, probs = 0.5),
            upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 1500x3, now 500x4]
#> mutate: new variable 'perc_change' (double) with 500 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped

Simulated annual condition factor by sub-division

div_index_list <- list()

for(i in unique(pred_grid$sub_div)){
  
  pred_grid_sub <- pred_grid %>% filter(sub_div == i & depth > 10 & depth < 110)
  
  pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 500)
  
  weight_sum_sub <- pred_grid_sub %>% 
    group_by(year) %>%
    summarise(density_cod = sum(density_cod))
  
  index_avg_sim <- get_index_sims(pred_avg_sim_sub,
                                  area = pred_grid_sub$density_cod) %>% 
    left_join(weight_sum_sub) %>% 
    mutate(est = est / density_cod,
           lwr = lwr / density_cod,
           upr = upr / density_cod) %>% 
    mutate(sub_div = i)
           
  div_index_list[[i]] <- index_avg_sim
  
}
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

div_index <- bind_rows(div_index_list)

Condition factor over time

# Spatiotemporally standardized indices
div_index2 <- bind_rows(div_index %>% dplyr::select(year, est, lwr, upr, sub_div) %>% mutate(sub_div = as.factor(sub_div)),
                        mutate(index_avg_sim, sub_div = as.factor("Weighted index")))

# Add in raw mean of data:
div_index2 <- bind_rows(div_index2,
                        d %>%
                          group_by(year) %>%
                          summarise(est = exp(mean(log_le_cren))) %>%
                          mutate(sub_div = "Empirical mean"))

# Add fixed effect predictions and the un-weighted index
div_index2 <- bind_rows(div_index2,
                        index_avg_sim_fe %>% mutate(sub_div = "Weighted index (fixed effects, year omitted)"),
                        index_avg_sim_uw %>% mutate(sub_div = "Un-weighted index",
                                                    upr = NA,
                                                    lwr = NA))


# Add sims for total and fixed effect prediction
index_comp_sims <- bind_rows(index_avg_sims_fe %>% mutate(sub_div = "Weighted index (fixed effects, year omitted)"),
                             index_avg_sims %>% mutate(sub_div = "Weighted index")) %>%
  filter(.iteration < 26)

# Plot!
p1 <- div_index2 %>% 
  filter(sub_div %in% c("Weighted index", "Un-weighted index", "Empirical mean", "Weighted index (fixed effects, year omitted)")) %>% 
  mutate(plot_facet = "Total and data") %>% 
  mutate(plot_group = ifelse(sub_div %in% c("Empirical mean", "Un-weighted index"), 1, 0)) %>% 
  mutate(line_group = ifelse(sub_div %in% c("Empirical mean", "Un-weighted index"), "1", "0")) %>%
  ggplot(aes(year, est, color = reorder(sub_div, plot_group))) +
  labs(y = "Le Cren's condition factor", x = "", color = "") +
  geom_ribbon(aes(x = year, y = est, ymin = lwr, ymax = upr, fill = sub_div), color = NA, alpha = 0.25) +
  geom_line(data = filter(index_comp_sims, .iteration < 26),
            aes(year, est, group = interaction(.iteration, sub_div), color = sub_div), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(aes(linetype = line_group), size = 0.6) +
  facet_wrap(~plot_facet, ncol = 1) +
  scale_color_brewer(palette = "Set2") +
  scale_fill_brewer(palette = "Set2") +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  theme_plot(base_size = 10) +
  theme(plot.margin = unit(c(0.4, 0.4, 0, 0.4), "cm"), 
        legend.position = c(0.24, 0.11),
        legend.direction = "horizontal",
        legend.spacing.y = unit(0, 'cm'),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(0.5, "cm"),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 8), 
        legend.background = element_rect(fill = NA)) +
  guides(fill = "none",
         linetype = "none",
         color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5,
                              override.aes = list(linetype = c(2, 2, 1, 1))))
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

p1
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf


p2 <- div_index2 %>% 
  filter(!sub_div %in% c("Weighted index", "Un-weighted index", "Empirical mean", "Weighted index (fixed effects, year omitted)")) %>% 
  mutate(plot_facet = "Subdivision") %>% 
  ggplot(aes(year, est, color = sub_div, fill = sub_div)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(y = "Le Cren's condition factor", x = "Year") +
  geom_line() +
  facet_wrap(~plot_facet, ncol = 1) +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2") +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  theme_plot(base_size = 10) +
  theme(plot.margin = unit(c(0, 0.4, 0.4, 0.4), "cm"), 
        legend.position = c(0.92, 0.88),
        legend.direction = "horizontal",
        legend.spacing.y = unit(0, 'cm'),
        legend.key.size = unit(0.8, "cm"),
        legend.key.width = unit(0.3, "cm"),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 8), 
        legend.background = element_rect(fill = NA)) +
  guides(color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5),
         fill = "none")

p2


p1 / p2 + plot_annotation(tag_levels = "A")
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf

#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf


# It gets messed up when I save in a knit, so i hashtag it here
#ggsave("figures/cond_index.pdf", width = 17, height = 20, units = "cm")

Coefficients

# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
                         filter(term %in% c("sigma_O", "sigma_E")),
                       
                       tidy(mfull, effects = "fixed", conf.int = TRUE)  %>% 
                         filter(!grepl('year', term))) %>%
  
  mutate(term = factor(term))

mfull_est <- mfull_est %>% 
  mutate(group = "Herring",
         group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
         group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
         group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
         group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
         group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
         group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
         group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
         group = ifelse(term == "sigma_E", "Random", group),
         group = ifelse(term == "sigma_O", "Random", group)
         )

mfull_est <- mfull_est %>% 
  mutate(scale = "Subdivision",
         scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
                                    "density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
         scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
                                    "depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
         scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
         scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
         )

# Quick calculation:
mfull_est %>% 
  mutate(par_type = ifelse(term %in% c("sigma_O", "sigma_E"), "re", "fe")) %>% 
  group_by(par_type) %>% 
  summarise(mean_est = mean(abs(estimate))) %>% 
  pivot_wider(names_from = par_type, values_from = mean_est) %>% 
  mutate(ratio = re/fe)


# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#>  [1] "biomass_her_sc"         "biomass_her_sd_sc"      "biomass_spr_sc"        
#>  [4] "biomass_spr_sd_sc"      "density_cod_rec_sc"     "density_cod_sc"        
#>  [7] "density_fle_rec_sc"     "density_fle_sc"         "density_saduria_rec_sc"
#> [10] "density_saduria_sc"     "depth_rec_sc"           "depth_sc"              
#> [13] "oxy_rec_sc"             "oxy_sc"                 "sigma_E"               
#> [16] "sigma_O"                "temp_rec_sc"            "temp_sc"

# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])

# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>% 
  mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))

ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_point(size = 2.5) +
  scale_color_manual(values = pal2, name = "") +
  scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
  scale_fill_manual(values = rep("white", 9), name = "") +
  labs(y = "Estimate", x = "Standardized coefficient") +
  scale_x_discrete(breaks = levels(mfull_est$term),
                   labels = c(expression(Herring[rec]),
                              expression(Herring[sub]),
                              expression(Sprat[rec]),
                              expression(Sprat[sub]),
                              expression(Cod[rec]),
                              expression(Cod[haul]),
                              expression(Flounder[rec]),
                              expression(Flounder[haul]),
                              expression(Saduria~entomon[rec]),
                              expression(Saduria~entomon[haul]),
                              expression(Depth[rec]),
                              expression(Depth[haul]),
                              expression(Oxygen[rec]),
                              expression(Oxygen[haul]),
                              expression(sigma[E]),
                              expression(sigma[O]),
                              expression(Temp[rec]),
                              expression(Temp[haul])
                              )) +
  coord_flip() +
  theme_plot() +
  theme(plot.margin = unit(c(0.2, 0.2, 0.4, 0.2), "cm")) +
  guides(color = "none", fill = "none", shape = guide_legend(ncol = 2))


ggsave("figures/effect_size.pdf", width = 17, height = 12, units = "cm")

# Just a test to see the labels were alright
ggplot(mfull_est, aes(reorder(term, estimate), estimate))

  
ggplot(mfull_est, aes(term, estimate)) +
  geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
  geom_point(size = 2) +
  labs(x = "", y = "Standardized coefficient") +
  coord_flip()


# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")

Predictions in space

pred_map <- predict(mfull, newdata = pred_grid)

pal <- wes_palette("Zissou1", 21, type = "continuous")

plot_map_raster +
  geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 8.5, dpi = 600)

# And a smaller map for selected years
lab_size <- 1.7

plot_map_raster +
  geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
              aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") + 
  facet_wrap(~ year, ncol = 2) +
  theme(legend.key.height = unit(0.1, "cm"),
        legend.key.width = unit(0.5, "cm"),
        legend.position = c(0.9, .036),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA),
        legend.text = element_text(size = 6),
        legend.title = element_text(size = 7.4)) +
  geom_sf(size = 0.1) + 
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)


ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")

Plot spatial and spatiotemporal random effects:

# Plot spatiotemporal random effect
plot_map_raster +
  geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 8.5, dpi = 600)

# Plot spatial random effect
plot_map_raster +
  geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
  scale_fill_gradient2() +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)

Calculate the “spatial trend” from the estimates:

# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
  mutate(id = paste(X, Y, sep = "_")) %>%
  split(.$id) %>%
  purrr::map(~lm(est ~ year, data = .x)) %>%
  purrr::map_df(broom::tidy, .id = 'id') %>%
  filter(term == 'year')

# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
  separate(id, c("X", "Y"), sep = "_") %>%
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

plot_map_raster +
  geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
  scale_fill_gradient2(midpoint = 0, name = "Slope (condition factor~year)") +
  geom_sf(size = 0.1)


ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)

Visualize conditional effects

Depth

# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))

nd_dep <- nd_dep %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)

cond_dep <- ggplot(p_cond_dep, aes(depth_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Depth[haul]))

Oxygen

# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))

nd_oxy <- nd_oxy %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         #oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)

cond_oxy <- ggplot(p_cond_oxy, aes(oxy_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se)))+
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Oxygen[haul]))

Temperature

# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))

nd_tem <- nd_tem %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         #temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)

cond_temp <- ggplot(p_cond_tem, aes(temp_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Temperature[haul]))

Cod

# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))

nd_cod <- nd_cod %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         #density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)

cond_cod <- ggplot(p_cond_cod, aes(density_cod_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Cod[haul]))

Sprat

# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))

nd_spr <- nd_spr %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         #biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)

cond_spr <- ggplot(p_cond_spr, aes(biomass_spr_sd_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.65, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Sprat[sub]))

Saduria

# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_rec_sc = seq(min(d$density_saduria_rec_sc),
                                                  max(d$density_saduria_rec_sc), length.out = 100))

nd_sad <- nd_sad %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         #density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model
p_cond_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)

cond_sad <- ggplot(p_cond_sad, aes(density_saduria_rec_sc, exp(est),
                                   ymin = exp(est - 1.96 * est_se),
                                   ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.75, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Saduria~entomon[rec]))

Plot together

# cond_dep + cond_oxy + cond_temp + cond_cod + cond_spr + cond_sad
#   plot_annotation(tag_levels = 'A') 
# 
# ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)

p_cond_dep2 <- p_cond_dep %>%
  mutate(variable = "Depth (haul)") %>% 
  dplyr::select(est, est_se, depth_sc, variable) %>% 
  rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_oxy2 <- p_cond_oxy %>%
  mutate(variable = "Oxygen (haul)") %>% 
  dplyr::select(est, est_se, oxy_sc, variable) %>% 
  rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_tem2 <- p_cond_tem %>%
  mutate(variable = "Temperature (haul)") %>% 
  dplyr::select(est, est_se, temp_sc, variable) %>% 
  rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_cod2 <- p_cond_cod %>%
  mutate(variable = "Cod (haul)") %>% 
  dplyr::select(est, est_se, density_cod_sc, variable) %>% 
  rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_spr2 <- p_cond_spr %>%
  mutate(variable = "Sprat (sub)") %>% 
  dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>% 
  rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_sad2 <- p_cond_sad %>%
  mutate(variable = "S. entomon (rec)") %>%
  dplyr::select(est, est_se, density_saduria_rec_sc, variable) %>% 
  rename(value = density_saduria_rec_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

conds <- bind_rows(p_cond_dep2, p_cond_oxy2, p_cond_tem2, 
                   p_cond_cod2, p_cond_spr2, p_cond_sad2)

ggplot(conds, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
                  ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  facet_wrap(~variable, scales = "free_x") +
  labs(y = "Le Cren's condition factor",
       x = "Scaled value") +
  theme_plot()


ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)

Condition model

Max Lindmark, Sean C. Andersson, Mayya Gogina and Michele Casini

18 May, 2023

Fit condition model with environmental and biological predictors

Load packages

Code for map plots

Read data

Read the prediction grids

Make spde mesh

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

Fit main model

Calculate approximate \(R^2\) manually using a modifed `r2.sdmTMB`

Plot residuals

Annual index using `get_index_sims`

Get sims from prediction with only fixed effects (-year) and random effects

Difference in Le Cren condition factor between the end and start of time period using sims

Simulated annual condition factor by sub-division

Condition factor over time

Coefficients

Predictions in space

Plot spatial and spatiotemporal random effects:

Calculate the “spatial trend” from the estimates:

Visualize conditional effects

Depth

Oxygen

Temperature

Cod

Sprat

Saduria

Plot together

Condition model

Max Lindmark, Sean C. Andersson, Mayya Gogina and Michele Casini

18 May, 2023

Fit condition model with environmental and biological predictors

Load packages

Code for map plots

Read data

Read the prediction grids

Make spde mesh

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

Fit main model

Calculate approximate \(R^2\) manually using a modifed r2.sdmTMB

Plot residuals

Annual index using get_index_sims

Get sims from prediction with only fixed effects (-year) and random effects

Difference in Le Cren condition factor between the end and start of time period using sims

Simulated annual condition factor by sub-division

Condition factor over time

Coefficients

Predictions in space

Plot spatial and spatiotemporal random effects:

Calculate the “spatial trend” from the estimates:

Visualize conditional effects

Depth

Oxygen

Temperature

Cod

Sprat

Saduria

Plot together

Calculate approximate \(R^2\) manually using a modifed `r2.sdmTMB`

Annual index using `get_index_sims`